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Metastable structures in macromolecular and colloidal systems are non-equilibrium states that often have 
long lifetimes and cause difficulties in simulating equilibrium. In order to escape from the long-lived metastable 
states, we propose a newly devised method, molecular Monte-Carlo simulation of systems connected to 3 reser- 
voirs: chemical potential p, pressure P, and temperature T. One of these reservoirs is adjusted for the ther- 
modynamic equilibrium condition according to Gibbs-Duhem equation, so that this adjusted 3rd reservoir does 
not thermodynamically affect phases and states. Additional degrees of freedom, i.e. system volume V and the 
number of particles N, reduce kinetic barriers of non-equilibrium states and facilitate quick equilibration. We 
show globally-anisotropic defect-free ordered structures, e.g. string-like colloidal assembly, are obtained via 
our method. 



Metastable states, e.g. amorphous solids of colloids Q 
and stalks in membrane fusion |2|, correspond to the regions 
where the free energy takes local minima in phase space. 
In conventional simulation of the canonical ensemble (NVT- 
ensemble), once the system is trapped in these metastable re- 
gions, the system tends to stay long over acceptable compu- 
tational time in these non-equilibrium states. Even if the sys- 
tem is allowed to leave non-equilibrium to ordered equilib- 
rium states, constant N and V cause defects in ordered struc- 
tures. This situation breaks global anisotropy of the ordered 
structure and compels the periodicity of the structure to fit the 
system size. For finding defect-free ordered structures, N and 
the rectangular system box size (l x , L y , L z ) need fine tuning to 
both the anisotropy and the periodicity, which are not known 
a priori. 

In advanced simulation techniques, e.g. multicanonical en- 
semble method [3|, the whole phase space is nearly homoge- 
neously sampled via artifitial weights that reduce occurrence 
probability of metastable states at constant N and \L X , Ly, L^j. 
However, the equilibrium microstates obtained by these ad- 
vanced techniques are limited to those with the given set of N 
and [L x , L, , LA. Free energy landscapes at the same particle 

density, N/V, and different sets of N and (l x ,L x ,L^ are not 
searched. These extensive variables need manual and simulta- 
neous fine tuning for the sake of determining the most stable 
state over these landscapes, which is a challenging task. For 
example, for perfect crystals, both N and \L X , L y , L,) must be 
integer multiples of the unit structure. In addition, prior to 
production runs, these techniques also require advanced pro- 
gramming and precise adjustments of the artifitial weights. 
Furthermore, unphysical sampling processes disallow tracing 
physical trajectories in the phase space. 

In the present article, we devise and perform molecular 
Monte Carlo simulation of systems connected to the 3 reser- 
voirs (3-reservoirs method), in order to study equilibrated 
colloidal assembly. We venture to connect the 3 reservoirs, 
though Gibbs-Duhem equation (GD eq.) limits the num- 
ber of reservoirs to 2. Additional degrees of freedom, N 
and \L x ,Ly, LA, are equivalent to the additional dimensions 



of the phase space and provide shortcuts, through the ex- 
tended regions of the phase space, from the non-equilibrium 
to the equilibrium state. The system itself is allowed to 
make the spontaneous and simultaneous fine tuning of N and 
\L X , L y , L.) for the ordered structure. Owing to physical sam- 
pling processes, the physical trajectories in the phase space 
are traced in simulation via our method. 

Guggenheim formally introduced Boltzmann factor of the 
ensemble with the 3 reservoirs [4|. Later, Prigogine and Hill 
also analytically studied this ensemble 0|6l. However, these 
early authors concentrated on formalism, i.e. mathematical 
aspects of the partition function, since their goal was to ob- 
tain a universal and generalized formula for a partition func- 
tion that is applicable to any thermodynamically acceptable 
ensembles [7 8 |. Physical aspects of this ensemble were to- 
tally left for the future. Here we study the physical aspects 
of this ensemble, mostly in intuitive and thought experimental 
ways. 

Despite GD eq., a system with the 3 reservoirs is readily 
constructed in experiments, for example, when a wall of a sys- 
tem box placed in the grand canonical ensemble is replaced 
with a free piston connected to another reservoir, i.e. P. We 
theoretically build thermodynamics and statistics of the sys- 
tems with the 3 reservoirs in the present article. We tentatively 
call this ensemble /iPF-ensemble. 

As an example of thermodynamics, we pack particles into 
a diathermal box with a free piston, placed in environment at 
constant (T,P). These two intensive variables determine the 
other intensive variables of this system, e.g. p, p — N/V, and 
free energy per particle. This means thermodynamic degrees 
of freedom of this system equal 2, which results from GD 
eq. In conventional simulation, N is also fixed at some value 
(A^Pr-ensemble), whereas states and phases of the system are 
independent of N; N only scales extensive variables. Briefly, 
phase diagrams constructed in Pr-plane are independent of 
extensive variables. Rather than fixing insignificant N, we 
connect a reservoir of p, determined from (T, P). This reser- 
voir does not affect phases and states, provided p and (T, P) 
are in equilibrium. This condition corresponds to a thermo- 
dynamically stable point (TSP) based on GD eq., resulting in 
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an equation of state linking T,P, and p. At TSP, extensive 
variables are freely scaled, i.e. indeterminate and fluctuating, 
while all the intensive variables are kept. This means, in sim- 
ulating this system, we can choose simulation runs at small N, 
which are computationally advantageous. Gibbs free energy 
per particle, which should be minimized in A^fT-ensemble 
prior to adding the 3rd reservoir of constant p., is still un- 
changed even after this 3rd reservoir is added. 

In the above example, (T, P) is given from the outside of the 
system; p is adjusted according to these (T, P) as an additional 
reservoir. Two other combinations (T,p) and P, and (p,P) 
and T also work. Gibbs free energy per particle, grand po- 
tential per volume, and the thermodynamic potential of pPS - 
ensemble per S , i.e. (E - pN + PV)/S , are simultaneously 
minimized in pfT-ensemble, where S and E denote entropy 
and internal energy respectively. This simultaneous minimiza- 
tion of the 3 free energy densities is originated from GD eq. 
When we connect the 3rd reservoir to the system and set the 
system at TSP, prior to the connection, the corresponding 3rd 
intensive variable needs adjusting. Through this adjustment, 
the corresponding free energy density is minimized. The same 
system at the same TSP is also built by the two other combina- 
tions of the 3 intensive variables, which results in the simulta- 
neous minimization of the two other free energy densities. On 
the other hand, in the other ensembles, any sets of correspond- 
ing 3 external parameters, e.g. (T,V,N) in A^Vr-ensemble, 
can be chosen without the adjustment. Therefore, only ther- 
modynamic potential of the ensemble, e.g. Helmholtz free 
energy F in A^V^-ensemble, is minimized. 

Here we sketch a system with the 3 reservoirs. Reservoirs 1 
and 2, composed of ideal particles with densities p£es an d PrlL 
specify p and P respectively. The system and these 2 reser- 
voirs are connected to reservoir 3, a thermostat. As another 
example of thermodynamics, we consider a system consisting 
of ideal particles. At TSP, a relation, p''s = Pies = P> holds. 
However, when p^s > Pres> both N and V diverge, since reser- 
voir 1 continues supplying particles and reservoir 2 supplying 
the volume. When p^ < p^s, both N and V vanish. There- 
fore, the system reaches equilibrium only at TSP. This also 
applies to systems of interacting particles. We utilize these di- 
vergence and vanishment as criteria for heuristically, e.g. by 
bisection method, determining TSP in our simulation. 

The detail of the statistical mechanics of 3-reservoirs 
method will be discussed in our forthcoming article (9). The 
partition function of pPr-ensemble is, 
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where kgT denotes thermal energy and Z^(T, V) the partition 
function of the system in A^VF-ensemble. The first pref actor 
of this equation, P/ksT, comes from conventions of NPT- 
ensemble ifTOll . Q.(T,P,p) is obtained in a similar manner as 
for the partition functions of the grand canonical and NPT- 



ensembles [10], as natural extension of these two ensembles. 
Boltzmann factor defined from D.(T,P,p) is, at fixed N, con- 
sistent with Boltzmann factor of iVPr-ensemble and, at fixed 
V, Boltzmann factor of the grand canonical ensemble. 

3-reservoirs simulation method is constructed based on 
conventional Monte Carlo simulation methods in the grand 
canonical ensemble (pVr-ensemble) and A^Pr-ensemble. 
Our algorithm is similar to Gibbs ensemble technique iflOl . 
which is used for simulating phase equilibria in NVT- 
ensemble. In one simulation step of 3-reservoirs method, the 
following processes are performed. 

i) with probability pc/2, trial particle insertion into the sys- 
tem 

ii) with probability pc/2, trial particle deletion from the sys- 
tem 

iii) with probability py, trial system size change 



iv) with probability 1 
particle 



Pg - Pv, trial displacement of one 



are chosen, where pc and py are constants fixed in intervals 
< pcPv < 1. 

Via conventional algorithms of trial particle insertion and 
deletion in pVr-ensemble ifTOl . we perform steps i) and ii), 
during which the system size is fixed. Steps i) and ii) sat- 
isfy the detailed balance condition |9|, since this condition is, 
due to the consistency between the Boltzmann factors of the 
present system and pVr-ensemble, consistent with the condi- 
tion in pVr-ensemble which has already been confirmed IflOl . 

Via conventional algorithms in iVPr-ensemble [10|, step 
iii) is performed, during which N is fixed. Step iii) also satis- 
fies detailed balance due to the consistency between the Boltz- 
mann factors of the present system and iVPr-ensemble |9)- 
However, unlike conventional Monte Carlo (MC) simulations 
in A^Pr-ensemble based on McDonald's method ITOl . L x ,L y , 
and L, are independently changed in our 3-reservoirs method. 

By Metropolis algorithm in A^VF-ensemble, step iv) runs. 

Here we discuss detailed balance and ergodicity, i.e. statis- 
tical nature, of 3-reservoirs method. Steps i) and ii) change N 
according to the detailed balance condition which corresponds 
to pVr-ensemble. Step iii) varies the system size accord- 
ing to the detailed balance condition corresponding to NPT- 
ensemble. The particle coordinates are updated in step iv), 
according to the standard Metropolis algorithm. As a result, 
3-reservoirs method fulfills detailed balance and ergodicity, 
based on the conventional ensembles which satisfy ergodicity. 
This means N, \L x ,Ly,LA, and the particle coordinates are, 
in the phase space, simultaneously searched and tuned to the 
equilibrium state. This also indicates that statistics of pPT- 
ensemble contradicts none of the three underlying ensembles 
with 2 reservoirs, which is consistent with the thermodynamic 
consideration. 

Maximization of statistical entropy per volume in pPT- 
ensemble corroborates the above results of statistics. The sta- 
tistical entropy is defined as —ks Yij Pj m Pj O |TT|, where the 
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suffix, j, denotes microstates of the system and pj is corre- 
sponding occurence probability. Due to indetermination of 
extensive variables in yufT-ensemble, we take volume density 
of extensive variables. The statistical entropy per volume is, 



Pjlnpj 



V; 



(2) 



The probability distribution, {p\,p2, . . . , Pj, . . .}, which max- 
imizes eq. is determined under constraints: 



z, 



(3) 



Tjj pj N v j) = < £ / y > • Tij p i N v j) = <*>>■ (4) 

where (• • • ) denotes the thermal average that should be spec- 
ified by the reservoirs. Equation (|3]l represents normalization 
condition. Equation (|4]) comes from the thermodynamic de- 
grees of freedom, which equal 2. Using Lagrange multipliers, 
we can solve this maximization problem 111 11 . The solution 
is0: 

pj = e- 1 exp [- (P/k B T) Vj - (l/k B T) Ej + ( M /k B T) Nj\ , (5) 
S/V = (1/T) \P + (E/V) - n (N/V) + k B T {\IV)\ . (6) 

Pj is equivalent to the partition function eq. ([T}. The thermo- 
dynamic potential of //PF-ensemble, 



<D = E - TS - fiN + PV, 



(7) 



is also obtained J9): = -k B T Ine, which is essentially zero 
compared with the other extensive variables in eq. |7]) in the 
thermodynamic limit. This result is consistent with the Eu- 
ler equation in thermodynamics. Moreover, the last term of 
eq. (|6} monotonically decreases with increasing V for finite 
V. This means that the principle of maximizing entropy re- 
stricts the system to the finite size for computer simulations 
that treat finite system size. These results are also confirmed 
by maximization problems of S /N and S IE. 

Assuming that the ensemble averages of extensive vari- 
ables, e.g. N and E, were determined, Guggenheim for- 
mally introduced Boltzmann factor (statistical weight) of 
//fT-ensemble, based on analogy between other conventional 
ensembles [4, 8|. However, Guggenheim's assumption cor- 
responds to, in the above calculation of the maximization 
problem, keeping the averages (E) and (N) fixed, instead of 
eq. (RJ. This contradicts the indetermination of extensive vari- 
ables, as was pointed out by Prigogine and Sack |5] |7). Pri- 
gogine showed that the resulting partition function diverges 
and therefore concluded that no physical meaning is found 
in this partition function J3] UJ . Free energy, O, statisti- 
cally determined via such partition function could be indefi- 
nite. Furthermore, in thermodynamics, <t> identically equals 
zero. Physical quantity that dominates juPr-ensemble, as cor- 
responds to F in A^Vr-ensemble, has been veiled since these 
early works. 

On the other hand, in the present work, eq. |5]) and cD un- 
der the constraints eq. (HI) indicate the partition function equal 



to e. This <D is smaller than the other extensive variables 
in eq. |7]) and vanishes in thermodynamic limit. Moreover, 
we have found, in discussion on thermodynamics, that NPT, 
fiVT, and p.PS -ensembles underlie yuPr-ensemble and that 
the 3 corresponding free energy densities of these underly- 
ing ensembles, rather than cD, are simultaneously minimized. 
This corresponds to the minimization of F in A^VF-ensemble. 

As an example of 3-reservoirs method, we simulate equi- 
librium states of model polymer-grafted colloids. The model 
colloids, 0"i in diameter, are interacting via spherically sym- 
metrical repulsive square-step potential with a rigid core |[T2"1 : 
0(r) = oo (r < tri), 4>(r) = e (o"i < r < cr 2 ), 0(r) = (<x 2 < 
r). The distance between centers of the particles is denoted 
by r and positive constants <r 2 and eo are the diameter and 
the height of the repulsive step resulting from the polymers 
grafted onto the colloidal particles. Simulating particles in- 
teracting via <p(r) in A^VF-ensemble, we have studied phase 
behavior of these colloidal systems lfl2l . These MC simu- 
lation results show that, at low T, high P, and cr%lcr\ ~ 2, 
our particles self-assemble into string-like assembly, an amor- 
phous solid. This assembly has also experimentally been ob- 
served [ 1 3 1 . 

In these recent studies at finite T in both 2 and 3- 
dimensions, the string-like assembly has been found to 
be metastable, locally directed in the same direction and 
globally-isotropic, in jVVr-ensemble. Equilibrium states at 
finite T are still open, though various ground states of our 
model at zero T have been found via genetic algorithms fl4l . 

Using 3-reservoirs method, at finite T, we simulate equili- 
brated states of our model system in 2-dimensions. cr\ and eo 
are taken as unit length and unit energy respectively. Dimen- 
sionless chemical potential is defined as ifTSl . 



li :=W^7 , -ln(A 2 / ( r?) = ln(p^ ( r?), A = 



\l2nmk B T 



where /; denotes Planck's constant, and m mass of the particle. 

In the initial state, iVo = 1254 particles are arranged on 
homogeneous triangular lattices in a square system box with 
the periodic boundary condition. Initially, pcr^ = 0.451. In 
step iv), a particle is given uniform random trial displacement 
within a square 0.40-; long in each direction. In step iii) |f9l , 
(L x ,L y ) is changed to ( L' x = L x + AU\ - 2£ v ), L' y = L y + 
AL(\ — 2g y )), where AL is constant length, fixed at O.Olcr], 
and % x and g y are random numbers uniformly distributed over 
intervals < £ t ,£ v < 1. p G = 10% and p v = 1 /No. With this 
pv, computational time is about twice as long as simulation 
in A^Vr-ensemble. In the following, we define 1 Monte Carlo 
step (MCS) as A^o simulation steps, o^/ch = 2.0 is fixed. 

First, simulation results at low temperature, k B T '/eo = 0.12, 
are discussed. In A^Vr-ensemble at this k B T/eo, the string- 
like assembly is observed in regions 0.412 ^ pcr\ % 0.481 
and the string length diverges at p(j\ « 0.451 [12]. Snap- 
shots of the system in the present simulation are shown in 
Figs. [TJa) and (b). Despite the different parameters, both 
the snapshots present globally-anisotropic defect-free string- 
like assembly, although a snapshot of the system simulated 
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in NVr-ensemble at pcr\ - 0.451 and N = 1200, shows, in 
Fig. [TJc) fl2l . globally-isotropic string-like assembly. Time 
evolution of p and N are plotted in Fig. [2] The same assem- 
bly is also observed at all the various (p' , Pcr^/eo) chosen in 
Fig. [2] Time evolution of p, Fig. |2ja), all fluctuates in the 
vicinity of pcr\ ~ 0.451, regardless of (p! ,Pcr\leo). These 
results indicate these parameter sets (p',P<r^/eo) are located 
in the vicinity of the same TSP, and that the system reaches 
the identical equilibrium state. Simulations started from dif- 
ferent initial conditions, e.g. different initial p and L x /L y , also 
reach these results. Simulations resumed from the instanta- 
neous state of Fig. [TJa) by disconnecting reservoirs 1 or 2, 
i.e. simulations resumed in NPT or yuVr-ensembles, keep the 
same assembly and p (data not shown), which corroborates 
our heuristic method for determining TSP. On the other hand, 
time evolution of N, Fig.[2|b), shows that the total system size 
depends on (p' ,Pcr^/ e ). N and V rise with large p' and small 
PcTj/e , whereas, regardless of the chosen (p',Pcr^/eo), com- 
putationally handy N and V are kept for a long time, within 
which good statistics of simulation results are obtained. With 
such results that are accurate enough to keep the equilibrated 
structure of the system in yuPr-ensemble, we are allowed to 
switch the ensemble to conventional ones, e.g. NPT or NVT- 
ensembles, and to perform long simulation runs, free from the 
divergence and the vanishment of N, of the equilibrated struc- 
ture obtained via the 3-reservoirs method. The system outside 
the chosen (p',Pcr^/e ), i.e. outside the vicinity of TSP, di- 
verges or vanishes quickly (data not shown), which facilitates 
heuristically determining TSP quickly. These results coincide 
with thermodynamics discussed in the introduction. Globally- 
anisotropic defect-free triangular crystals of the outer cores, 
0"2, are also observed at this temperature (data not shown) (9), 
which shows our method is applicable to crystalline states. 
Since any crystals fit into rectangular system boxes with peri- 
odic boundary conditions, other shapes of system boxes, i.e. 
a combination of Parrinello-Rahman technique IfTUl and our 
method, are left for the future. 

At k B T/e Q = 1.5, Po^/eo = 10 and \£ = 12.50, a disordered 
state is observed and at keT/eo = 0.01, Pcr^/eo = 1 and p.' = 
34.53, a triangular crystal (data not shown). 

In conclusion, we have shown that yuPr-ensemble is built 
by combining NPT, pVT, and pPS -ensembles. The 3 corre- 
sponding free energy densities are simultaneously minimized 
in //fT-ensemble, rather than <t> is. Unlike the early works 
shown above, thermodynamics outside TSP and maximiza- 
tion of S /V have also been discussed. We have devised 
3-reservoirs method and simulated the colloidal suspension 
which are interacting via hard core with step repulsive po- 
tential. Unlike other advanced techniques, our method allows 



tracing the physical trajectories in the phase space and quickly 
starting production runs without advanced programming and 
large amounts of preliminary simulation. This facilitates and 
reduces the total work flow of simulation studies, including 
data analysis. Simulation results coincide with thermodynam- 
ics. Simultaneously tuning N and the system size, which 
is feasible only via 3-reservoirs method, and crossing non- 
equilibrium, the systems reach equilibrium at low T . These 
results illustrate our method is applicable to ordered equilib- 
rium states of various physical systems at finite T. This unique 
advantage, i.e. the simultaneous tuning of N and the system 
size, of 3-reservoirs method could overturn previous simula- 
tion results obtained via the other simulation techniques. The 
equation of state linking p, P, and T is also calculated. 
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FIG. 1. Snapshots of the system, (a), k B T/e = 0.12, Pcr\le = 1 and /i' = 27.63. (b), k B T/e = 0.12, Pcr]/e = 1.1 and n' = 29.93. Black dots 
represent the centers of the particles and grey lines denote networks of overlaps between the particles. A snapshot of the system simulated in 
NVT-ensemble at Ar s 77e = 0. 12, per] = 0.451 and N = 1200 is also shown in (c) fl2l . 



(a) 




MCS/10' 



FIG. 2. Time evolution of p, (a), and time evolution of N, (b), sam- 
pled at k B T/e = 0.12 and various Pa\le , \i' . Lines 1) to 5) are 
results at Por]le = 1.0. 1): y.' = 27.63. 2): fj? = 29.93. 3): 
H' = 31.54. 4): // = 32.23. 5): n' = 34.53. Lines 6) to 8) are 
results at // = 29.93. 6): Pa]le Q = 0.95. 7): Pcr\le = 0.98. 8): 
Po- 2 Jeo = 1.1. 



